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<N i ABSTRACT 

The effect of magnetic shear and shear flow on local gravitationally induced 
\ instabilities is investigated. A simple model is constructed allowing for an arbi- 

m ' 

(N 

A transformation to shearing magnetic coordinates achieves a model with plasma 



flow along the magnetic field lines where the coordinate lines are coincident with 
the field lines. The solution for the normal modes of the system depends on two 



Q-V parameters: the Alfven Mach number of the plasma flow and the entropy gradi- 

ent. The behavior of the unstable normal modes of this system is summarized by 
a stability diagram. Important characteristics of this stability diagram are the 
following: magnetic shear is stabilizing and the entropy gradient must exceed a 
threshold value for unstable mode growth to occur; flow acts to suppress mode 
growth in a substantially unstable regime as expected, yet near marginal stabil- 
ity it can lessen the stabilizing effect of magnetic shear and enhance the growth 
rates of the instability; and, as the Alfven Mach number approaches one, the 
instability is completely stabilized. Analytical work is presented supporting the 
characteristics of the stability diagram and illuminating the physical mechanisms 
controlling the behavior of the model. A derivation of the stability criterion for 
the case without shear flow, asymptotic solutions in the limit that the Alfven 
Mach number approaches one and in the limit of zero growth rate, a complete 
WKB solution for large growth rates, an exactly soluble bounded straight field 
case, and energy conservation relations are all presented. The implications of 
this work for astrophysical and fusion applications and the potential for future 
research extending the results to include compressibility are discussed. 



Subject headings: MHD — gravitational instability — magnetic shear — velocity shear 
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1. Introduction 

The gravitational stability of a fluid against convective motion has been extensively 
studied over the past century. Pioneering examinations of the stability of unmagnetized 
and magnetized compressible fluids have been conducted by Schwarzschild (1906), Newcomb 
(1961), and Parker (1979). These papers have had a profound influence on diverse subjects 
from the dynamics of astrophysical objects to the confinement of plasma in a fusion device. 
Here we examine the effects of shear in the magnetic field and of an applied shear plasma flow 
on stability against gravitational interchange. These effects change the stability properties 
and our results are important for many applications. 

We construct a simple model to study the effects of magnetic shear and shear flow on the 
stability properties of a magnetized plasma in a gravitational field. We derive the equations 
which determine the behavior of this model in the Boussinesq limit. The equations depend 
on two parameters, the plasma flow Alfven Mach number and the entropy gradient. We 
conduct a numerical study of the normal modes of instability and summarize the behavior 
of the unstable modes by a stability diagram. The stability diagram demonstrates three 
important characteristics. First, the entropy gradient must exceed a threshold value for 
unstable mode growth to occur. Therefore, the shear magnetic field can stabilize a nonzero 
entropy gradient. Second, as expected, shear flow does act to suppress unstable mode growth 
when the system is at a substantially unstable point in parameter space. But, surprisingly, 
near marginal stability, shear flow actually enhances the growth rates of the instability and 
also lowers the threshold entropy gradient necessary for instability. The effect of magnetic 
shear — to stabilize the plasma and increase the threshold entropy gradient — is diminished by 
the addition of shear flow. The system can extract energy from the shear flow to further drive 
the system to instability. Third, as the Alfven Mach number approaches one, the unstable 
growth rate is suppressed; the normal modes of the instability are completely stabilized when 
the plasma flow exceeds the Alfven speed. Here, the unstable region in space where a mode 
can grow moves faster than any perturbation in the system; any disturbance will be swept 
downstream out of the unstable region, leaving behind a stable plasma. 

Analytical work helps us to understand the mechanisms responsible for these characteris- 
tics. First, the stability criterion for the case without shear flow is derived demonstrating that 
a shear magnetic field can support a positive entropy gradient. Next, asymptotic solutions 
demonstrate analytically that stabilization occurs as the Alfven Mach number approaches 
one, yet the threshold entropy gradient for instability goes to zero in the same limit. In 
addition, a complete WKB solution in the limit of a large growth rate demonstrates both 
the stabilization by flow at large growth rates and destabilization near marginal stability. 
Then, a bounded straight field case is solved exactly to show that the lowering of the thresh- 
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old entropy gradient with increased shear flow is a characteristic of plasma flow along the 
field lines and not dependent on the magnetic shear in the general model. Finally, energy 
conservation relations are derived and analyzed. 

In section 2, we describe the model under consideration and derive the governing system 
of equations. The numerical stability analysis for the general model and the stability diagram 
are presented in section 3.1. Section 3.2 contains the analytical results illuminating the 
characteristics of the stability diagram. We describe an exactly soluble, bounded, straight- 
field case in section 3.3. Finally, in section 4, the implications of this work on galactic physics, 
accretion disk physics, solar physics, and tokamak confinement are discussed. 



2. Derivation of Equations 

In this section, we derive the equations for linear perturbations of a vertically stratified 
atmosphere with shear magnetic field and shear flow in the high-/?, or Boussinesq, limit. We 
motivate and apply a coordinate transformation which casts the problem in its most natural 
form. Investigating the limit of the most unstable modes, we derive a system of three 
coupled first-order ordinary differential equations which capture the lowest order behavior 
of the model. 



2.1. Setup and Coordinate Transformation 

Consider a stationary state for an ideal plasma with mass density p(x) and thermal 
pressure p(x) and an embedded horizontal shear magnetic field given by 

B = B (z + ^-y) . (1) 

We then impose a shear flow on this plasma given by 

x 

v = v -y, (2) 
i v 

and include gravitational acceleration in the vertical direction given by g = — gSc. Equilib- 
rium force balance yields 

Figure 1 shows the geometry of this shear magnetic field as well as the imposed shear flow 
on the system. 
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The instabilities of this plasma are expected to have a short perpendicular wavelength 
and a long parallel wavelength (with respect to the magnetic field) so as to maximize motion 
in the vertical direction and minimize field line bending (Newcomb 1961). Short perpen- 
dicular wavelengths, however, are rapidly sheared apart by the perpendicular shear flow. 
We would thus like to transform to a coordinate system with two properties: first, that 
the flow is along the magnetic field lines; and, second, that field lines are coordinate lines. 
The application of sheared coordinate systems to simplify a problem of this nature is well 
documented. Roberts and Taylor (1965) employed a coordinate system in which the field 
lines are coincident with the coordinate lines to describe the Rayleigh- Taylor instability of a 
fluid supported by a shear magnetic field; Goldreich and Lynden-Bell (1965) used a shearing 
coordinate system to attack the problem of local gravitational instabilities in a system with 
shear flow arising from differential rotation. 

We transform the shear velocity to a parallel velocity by transforming to a frame moving 

in z 



i.e. parallel to B . We construct field line coordinates by transforming the y coordinate to 



The x coordinate is left unchanged, x' = x. It is easy to verify that x' and y' are constant 
along field lines (B • Vx' = B • Vy' = 0) and that B = B Vx' x Vy' . The surfaces of 
constant y' twist from vertical at z' — to almost horizontal as z' — > ±oo. This geometry is 
shown in Figure 1. 

At first, introducing the field line coordinates seems unhelpful since it introduces explicit 
z' dependence into the equations. The problem also has x dependence that arises from the 
variation of B, p(x), and p{x). Indeed, one way to tackle this problem is to Fourier transform 
in y and z and solve for the x dependence. However, the lowest order solution in the twisting 
coordinate system is a superposition of these Fourier modes, or, complementarily, a Fourier 
solution can be constructed by a superposition of these twisting modes (Roberts and Taylor 
1965). We summarize the relationship between these representations in Appendix A. We 
consider the solutions in the twisting coordinates to be more physically relevant since they 
are localized in z' . 



Z = Z + Vft. 

where Vf — l -f-v . In this transformed frame, the velocity becomes 



(4) 




y = y-~r- 



xz' 



(5) 
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2.2. Application of Ideal MHD 

The basic equations of ideal MHD include the momentum equation written in terms of 
the gradient of total pressure (thermal and magnetic), the magnetic tension force, and the 
gravitational force, 

Dv / B 2 \ B • VB 



the induction equation in the limit of zero resistivity, 

— = Vx(vxB), (7) 

the continuity equation, 

g = -,V-v. (8) 

and the adiabatic equation of state, 

Here = + v • V denotes the Lagrangian derivative, v represents the plasma velocity, 
and T is the adiabatic index. These equations must be evolved subject to the constraint that 

V-B = 0. (10) 

Equations (6)-(9), linearized about the equilibrium (equations (l)-(3)), yield: 

3 \ r . „ l n / t B 5Bn\ B d5B 5B VB g5p± M ^ 
i + v f —)5v + 5-vVv = VUp+— + - - +— o-^i- 11 



i + ] SB + 5v • VB = B -^ + <5B • Vv - B V • 5v 

( 7 ' + U/ ^ 7 ) 5p = ' Vp ° ~~ PoV ' ^ V 
( 7 ' + U/ ^ 7 ) Sp = ~ 6v ' Vp ° ~ rPoV ' 5v ' 



(12) 
(13) 
(14) 



where we have taken all quantities to vary in time as e 



In a straight field (Newcomb 1961) and a shear field (Roberts and Taylor 1965) without 
flow, the most unstable perturbations are incompressible to lowest order and have a small 
horizontal perpendicular wavelength. Such perturbations maximize vertical motion which 
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extracts energy from the gravitational potential energy and minimize horizontal motions 
which extract no energy. To isolate these motions, we impose the ordering 

£ = * ~ °(9 (15a) 
w ~ °W) (15b) 

(l5c) 

f ~ ^(e 1 / 2 ). (15d) 

where e = (his) -1 <d 1 is the ordering parameter or our problem. Clearly, all perturbed 
quantities vary as e lky ' . It is also convenient to define the vectors, 

b = % (16a) 

e A = Vy' (16b) 
W x Bp 

e± = b • ( 16c ) 

-DO 

The perturbed plasma velocity and magnetic field are projected along these directions, i.e. 

5v = Sv±e± + Sv A e A + (foyb 

5B = 5B ± e ± + 5B/\G/\ + 5B {{ h. 

Note that the basis vectors ej_, e A , and b are neither unit vectors nor constant in space — e.g., 
B ■ Ve A = — ^x. We expand all perturbed quantities in powers of e 1 / 2 and denote order as 

a superscript — for example, 5v± = Y^=o &v^e n l 2 . The ordered, perturbed quantities and 
operators are substituted into equations (11)— (14). Equations (12)-(14) at 0(e~ 1 ) and the 
e A projection of equation (12) yield 

Sv® = SB A 0) = 0. (17) 

Thus, the dominant motion is along the constant y' surfaces in the e_i_ direction. Equa- 
tions (12)-(14) at Oie' 1 / 2 ) produce 

V • (5v i ° ) e ± + Sv A 1] e^ = 0. (18) 

Thus, the perpendicular motion is incompressible to lowest order. At 0(e~ l ), the e A com- 
ponent of equation (11) gives 

5p (o) + II = o. (19) 

47T 
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Equation (19) expresses the fact that, on the time scales of interest, pressure balance is 
achieved across the convective eddies (in the e A direction). In a subsidiary expansion, we 
take the high-/?, or Boussinesq, limit {(3 = -pr)- Thus, equation (19) reduces to Sp^ = 0, 
and we find from equation (14) that 

5v (0) • Vp = -I>o (V • <Sv) (0) . (20) 

B 5B (0> 

At C(e -1 / 2 ), the e A component of equation (11) yields Sp^-\ = 0, or taking the high-/? 

limit, 5p^ = 0. The final stability equations are obtained from the sum and difference of the 
projections of equations (11) and (12) at 0(1) and from equation (13) using equation (20) 
to substitute for V • 5v. After some tedious but straightforward algebra, we obtain 

_ (1 _ m _± = - lA+ + (1 + M)—^ - JrT - w (2!) 

d A 7 <? 

(1 + M)-^ = -.A. - (1 - M)—A + + ITT - ¥7 , (22) 

, r d-S 1 -/i-U - 1 , - . 

M H = ~~< s - 2 (TTPF 54 ' (23) 

where A + and A_ are the Elsasser variables defined by 

6Bi + 6 ^] (24) 



+ (1 + ^ 2 )V2 \ Bo v A 



1 + z 1 ) 1 ' 2 V Bo 



and the entropy is given by 

s = 

U A VPO, 

We have normalized so that z = z'/l B and 7 = ^'Ib/va- Equations (21)-(23) contain 
two parameters: the Mach number of the plasma flow with respect to the Alfven speed 

( Va = (4^72)' 

M — — , (27) 

^A 



4 = ?(--^). (28) 



and the entropy gradient, 

gj| f Po go 
^a VPo Tp 0/ 

where the primes denote differentiation by x. With the boundary conditions that A + — > 0, 
A_ — > 0, and s — > as | z |— > ±00, equations (21)-(23) define an eigenvalue problem for 7. 
Solution of the stability equations yields 7(M, s' ). 
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Equations (21)-(23) have a simple physical interpretation. A+, the Alfven wave going in 
the negative z (upstream) direction, travels at the (normalized) speed 1 — M. The A + wave 
is coupled to the A_ wave by magnetic and velocity shear (the A_ term in equation [21]). 
The A + wave is driven by gravity via the s term in equation (21). A_, the Alfven wave going 
in the positive z (downstream) direction, travels at speed 1 + M, is coupled to A + , and is 
driven by s. The variable s is proportional to the density or entropy perturbation and it is 
driven by both Alfven waves, as seen in equation (23). 



3. Stability Analysis 

In this section, we discuss the unstable eigenvalues (Re 7 > 0) and eigenfunctions 
of equations (21)-(23). We have not examined the stable part of the spectrum in detail 
although numerical results indicate a continuum along the imaginary 7 axis. Two properties 
of equations (21)-(23) show that it is sufficient to examine stability in the region < M < 1. 
First, note that 7(— M, s' ) = 7(M, s' ) since we can map equations (21)-(23) onto themselves 
by the changes M — > — M , A + — > A-, A- — > A + , z — > —z, and s — > s. Second, note that 
the three asymptotic solutions as \z\ — > 00 are: 

72 / 1 7 2 \ 

A + ~ e~ s, A_ ~ eU-ei=Mj (29) 

A_ ~ e i+M s,A + ~ O (-ew j (30) 



s~e^r A + ,A_~ O(-e^r). (31) 



If M > 1, there are no acceptable (decaying) asymptotic solutions as z — > —00 for Re 7 > 0. 
Thus, M > 1 has no unstable eigenmodes. Physically, this is because all solutions, even the 
upstream propagating Alfven wave A + , are swept downstream. 

Without flow, MHD stability equations are self-adjoint (Bernstein et. al. 1958) and 7 2 
is real. With flow, no such property is known and 7 2 can be complex. However, in all our 
solution methods, 7 2 has been found to be real for this problem, although we have not been 
able to prove that this is rigorously true. The discrete positive real values for 7 correspond to 
unstable growing modes, and the continuum along the imaginary 7 axis represents traveling 
Alfven waves. 

In the following subsection, we present the numerical solution of equations (21)-(23). 
Various analytical limits that illuminate the numerical solutions are treated in subsection 3.2. 
An exactly soluble model with a straight magnetic field that demonstrates qualitatively 
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similar behavior is presented in subsection 3.3. Finally, in subsection 3.4, energy constraints 
on the instability are discussed. 

3.1. Numerical Solutions 

We looked for normal mode growth in the system defined by equations (21)-(23) using 
three different numerical methods. We directly solved for the eigenvalues 7 of this system 
by matrix solution of the corresponding finite difference equations using the commercial 
numerical routine package LAPACK. We also found the eigenvalues of the equations to high 
precision using a 1-D shooting code in z driven by an adaptive step-size, fourth order Runge- 
Kutta method with fifth order correction (RK45). Finally, for equations (21)-(23) with 7 
replaced by d/dt, an initial- value code employing Barton's method (Centrella and Wilson 
1984) for second order accuracy in time was written to determine the fastest growing mode 
for any choice of parameters. Results from all three codes were consistent. 

A stability diagram of unstable normal-mode growth rates and stable regions over the 
parameter space defined by M and s' Q is presented in Figure 2. As we have already noted, 
no growing mode exists for M > 1. As well, it is obvious that for the non-positive values of 
the entropy gradient, s' Q < 0, there can be no unstable mode growth since the atmosphere 
is stably or neutrally stratified; this is demonstrated by the energy arguments presented in 
section 3.4. Hence, Figure 2 need only cover the region of (M, s' ) parameter space defined 
by < M < 1 and s' > to include all possible unstable mode growth. 

Several features of Figure 2 are important to emphasize. First, for a system without 
plasma flow (along the line M = 0), we see that the entropy gradient must exceed a threshold 
value, s' > 1/4, in order to become unstable when the fluid is supported by a shear magnetic 
field. Second, the qualitative effect of increasing the plasma flow (increasing M) on the 
instability growth rate depends on both the growth rate and the plasma flow. Away from 
marginal stability (7 > 1), an increase of the plasma flow — equivalent to moving along a 
horizontal line to the right on the stability diagram — decreases the instability growth rate 
as one may expect. Here, the growing perturbation is sheared out horizontally so that less 
gravitational potential energy is extracted by motions along the constant y' surface, reducing 
the instability growth rate. But near marginal stability (7 < 1), an increase in the shear flow 
effects the instability growth rate differently depending on the relative orders of the growth 
rate, 7, and of one minus the Alfven Mach number, 1 — M. The diagonal dotted line in 
Figure 2 denotes 7 = 1 — M. In region (II) of the diagram, 7 > 1 — M and flow suppresses the 
instability. But, in region (I), where 7 < 1 — M, an increase in the flow actually enhances the 
instability growth. This unexpected result can be explained with some physical insight. The 
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point where the instability can grow is localized at z' — in our transformed coordinates; 
this is where the constant y' surface is vertical and motions along that surface can extract 
the most gravitational potential energy with which to drive the instability. When a shear 
plasma flow is introduced into the system, this is manifested in our transformed system by 
a plasma flow along the field lines in the +z' direction. This causes the Alfven modes in 
the +z' and —z' directions to propagate at different speeds in our transformed system. The 
counter-propagating mode (A + in our model) is slowed down, spending more time in the 
region around z' — where unstable mode growth occurs. Hence, the instability is enhanced 
by a shear flow in the plasma. The final point to be gleaned from Figure 2 concerns the 
behavior as M — > 1, or as the Alfven Mach number approaches one. In this region, every 
contour corresponding to a finite growth rate asymptotes to s' — > oo; hence, the instability 
is stabilized as the Alfven Mach number approaches one. 



3.2. Analytical Limits 

The stability diagram, Figure 2, prominently displays the three main characteristics dis- 
cussed in section 3.1: a threshold entropy gradient, s' > 1/4, necessary for instability in the 
system without plasma flow; for increasing plasma flow, a reduction of the unstable growth 
away from marginal stability, but an enhancement of that unstable mode near marginal sta- 
bility including a decrease in the threshold entropy gradient necessary for instability; and the 
stabilization of unstable normal modes as the Alfven Mach number approaches one. Each 
of these characteristics is relevant in a different region of the (M, s' Q ) parameter space of 
the diagram. By examining the model in each of these regions of parameter space, we can 
confirm and explain our results analytically. 

First, we examine the criterion for instability in the absence of plasma flow in sec- 
tion 3.2.1; this corresponds to the left vertical axis of the stability diagram where M — 0. 
Next, in section 3.2.2, we conduct an asymptotic analysis in the M — ► 1 limit — region (II) 
of the stability diagram — to show that the plasma is indeed stabilized as the Alfvenic Mach 
number approaches one. In section 3.2.3, we investigate the reduction of the threshold for 
instability with plasma flow; this covers the lower, right-hand side of region (I) of the stabil- 
ity diagram. A WKB analysis for a large instability growth rate, presented in section 3.2.4, 
yields the behavior of the system in the central and upper portion of the stability diagram; 
the suppression of the growth rate by flow in region (II) and its enhancement in region (I) 
are verified by the eigenvalue condition j(M, s' ) obtained in this analysis. 



- 11 - 



(1 + 2 >£ 



[ 7 2 (1 + z 2 ) - 4] v = (32) 



3.2.1. Stability Criterion without Flow 

Here we obtain the stability criterion for a magnetized fluid supported by a shear mag- 
netic field in the Boussinesq limit with no shear flow. By using the substitution v = 
Sv±/(1 + z 2 ), the equations without plasma flow (M = 0) can be simplified to a Sturm- 
Liouville equation of the form 

d_ 

dz 

over the interval (— oo, +oo). The boundary conditions on this system necessitate that v — > 
as z — > ±oo. From Sturm's First Comparison Theorem (Ince 1926), we know that, as the 
eigenvalue 7 2 is increased, the solution will oscillate less rapidly with zeros of the function 
v leaving the interval — oo < z < oo at the boundaries. Thus, if the solution with 7 2 = 
oscillates, we can increase 7 2 until the boundary conditions are satisfied, so there will be an 
unstable solution. Note also that the fastest growing mode has no zeros in the interval and 
must be even in z. 

Let w(z) satisfy equation (32) with 7 2 = 0. Substituting a series solution of the form 

oo 

Y^a n (l+z 2 )- (n+a) - (33) 



oo 

w(z) 

n=0 



(34) 



in equation (32) (with 7 2 = 0), we obtain the recurrence relation, 

a n _ 4(n + a- l) 2 

a„_i ~ 4(n + a - 1/4) 2 + s' - 1/4' 

with 

« = i ± i V / l^H- (35) 

When s' > 1/4, solution (33) oscillates for z — > oo, and clearly an unstable solution 
exists. Let us therefore consider stability when s' < 1/4. We take the positive sign in 
equation (35); then w 2 is integrable for z — > oo. We note that all a n are positive if a>o 
is positive. We choose a > such that w > 0. The series solution for w(z) given by 
equation (33) is non-differentiable at z = (as we see below). Thus, we cannot use the 
solution from equation (33) over the whole interval and must restrict its use to z > 0. Let us 
suppose (for contradiction) that there exists at least one unstable solution of equation (32). 
Further, let v (z) be the most unstable solution — as noted above, v (z) must be even in z 
and have no zeros in the interval — oo < z < oo. We therefore choose v (z) > everywhere. 
It is straight forward to show that 

7o / wv (l + z )dz, (36) 



o+ 
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where the limit + is infinitesimally above z — 0. Since the integral and t> in equation (36) 

dw I 
ciz lo+ 



are positive, we have stability, 7q < (a contradiction), if ^| n+ < 0. Since every term in 



the series equation (33) is a monotonically decreasing function of z we expect ^| 0+ < 0, 
but since the limit is nonuniform we take a more careful approach. We determine the sign of 
the limit ^| 0+ from an examination of the convergence of the series. It is straightforward 
to show that 

as n — > oo with A a positive constant. Thus, the series for w(z) (equation [33]) converges 
for z > 0. However, the series for ^ converges for z > but not for z — 0. Let us write for 

rfw _y^ -2^(n + a)a ra _ ^ A^n~ 1/2 

^ Z^(l + z 2)(n+a+l) - UZ Z Z^(1 + Z 2)™' ^ 



n=0 v 7 n=iV 

where C is a positive constant and N is a large number in the range 1 C <C z~ 2 . Using 
(1 + z 2 y n ~ e~" z2 , we obtain 



^ ~ -Cz - 2A J^Zn n-^ze-™ 2 ~ -Cz - 2A /~ n" 1 / 2 ^"^^ = -Cz - 2A f^ g e^dp 

~ -Cz - 2A J °° e^p = -Cz - AypR. (39) 

Thus, the limit of ^ as z — > + is — y/wA, ^.e. negative. From equation (36) we conclude 
that 7q < and there are no unstable modes for s' < 1/4. Thus, the necessary and sufficient 
condition for instability is Sq > 1/4. 

This criterion can also be written 

^ - A > I A (40) 

This confirms the result in Figure 2 — that a threshold value of the entropy gradient, given 
by Equation (40), must be exceeded in order to cause instability when the fluid is supported 
by a shear magnetic field. Clearly, without magnetic shear (Ib — > oo), the usual criterion, 
s' > 0, holds and the motion is the simple interchange of field lines. With magnetic shear, 
the field lines must be bent since the interchange of field lines is impossible with finite 
displacements — thus magnetic shear is stabilizing. 



3.2.2. Asymptotic Solution in the M — > 1 Limit 



An asymptotic, boundary layer analysis can be carried out in the limit that M — > 1. 
This asymptotic expansion is described in Appendix B. The eigenvalue condition derived 
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for this limit in that appendix is 

7 2 ~ s' (l - M 2 ). 

This relation explains the upturn towards infinity of the constant growth rate contours in 
region (II) of Figure 2 as M — > 1. 



3.2.3. Solution in the 7 — > 0, M — ► 1 Lzmii 

We demonstrate that the threshold entropy gradient for instability (the limit that 7 — > 
0) decreases as the plasma flow velocity is increased and that the threshold value of s' 
approaches zero linearly as M — > 1. Letting 1 — M ~ e, we can redefine the following 
variables in terms of e: 7 = 67, A_ = eA_, s = es and s' = es' Q . For M — > 1, we can 
drop terms of order e 2 and cancel e from each remaining term to yield the simplified set of 
equations 



= + -A ■ — - (41a) 



^+ + TT-Z^ ( 41b ) 



t = -kjiT^W 1 - (41c) 

These equations are now independent of e, or, equivalently, are independent of M. The 
equations will hold true for constant values of 7 and s' . In this case, if we have a negligibly 
small value 7 — > 0, we obtain the corresponding threshold value of s' = s' (l — M), where Sq 
is a constant. Thus, the threshold entropy gradient for stability must linearly approach zero 
as M — > 1 , as seen in the lower right-hand corner of Figure 2. 



3.2.4. Wiirtf Analysis 



A complete solution of the model for a large growth rate 7 can be constructed if we 
assume an ordering, for a small parameter e, such that 7 ~ (^(e^ 1 ), ^ ~ 0(e~ 1 ), s ~ (9(e _1 ), 
and s' ~ (9(e -2 ). In this case, we can neglect the second term on the right hand side of both 
equations (21) and (22). Converting back from Elsasser variables to 5B± and Sv± notation, 
combining the three equations into a single second order equation, and neglecting the term 
Msj- { jy^yj^ j (which can be shown to be small), yields 



(1-M 2 )-— ^-2 7 M 



dz 2 



dz 



7 2 " 



1 + z 2 



Svi 



0. 



(42) 
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Changing variables with an integrating factor to v — exp ( x-m' 1 ) ^ v± i we obtain 

1 



dz 2 



;i-m 2 ) 2 



7 2 - 



l + z 2 



v = 0. 



(43) 



We assume a WKB solution of the form e l f k ( z ) dz an( j find 



f = f exp 



1 - M 2 



s' (l-M 2 ) 



1 + r 



/2 



7 



1/2 



(44) 



The turning points are at z — ±z where 



^_ s' (l-M 2 ) 

T 



I. 



(45) 



For | z |> zo, the WKB solutions are decaying exponentials; in the region — z < z < z , the 
WKB solution is an oscillatory function. 

In the usual way (Bender and Orszag 1978), we obtain the Bohr-Sommerfeld quantiza- 
tion condition 



s'o(l ~ M 2 ) 
l + z 2 



nl/2 



7 



dz = 2nvr(l - M 2 ). 



The growth rate of the nth mode is then given, for small z , by 

7 2 = [s' - 2n(s' ) 1/2 (l - M 2 ) 1 / 2 ] (1 - M 2 ). 



(46) 



(47) 



For a large growth rate 7 ~ 0(e x ) and the lowest, nontrivial eigenmode n — 1, we can solve 
this eigenvalue condition for s' Q to obtain 



r 



1-M 2 



+ 27. 



(48) 



This condition agrees with the behavior of the constant growth rate contours in region (II) 
of Figure 2. Although we do not expect the eigenvalue condition in the WKB approximation 
to be precise in the limit of small growth rate, 7 ~ O(e), we do find to lowest order the 
qualitatively correct form, 

s' ~ 1 — M 2 , (49) 
that the contours in region (I) decrease like 1 — M 2 . 
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3.3. Bounded, Straight-Field Case 



The general model denned by this paper has the characteristic that the shear in the 
magnetic field localizes the region of instability around where the constant y' surface is 
vertical (this corresponds to z' — in our transformed coordinate system). With this 
characteristic as our guide, an exactly soluble, simplified model can be constructed which 
demonstrates the same qualitative behavior displayed in Figure 2. We construct a case with 
a straight magnetic field and a plasma flow along the field lines which has boundaries at 
z = ±L. Since the explicit z dependence drops out of the equations in the straight-field 
limit (Ib — > oo), we Fourier transform in z to obtain the algebraic dispersion relation 



If we remove the plasma flow from the system by setting M — 0, this dispersion relation 
agrees with the results of Newcomb (1961) for a straight-field without flow in the Boussinesq 
limit. 

We solve for the three solutions for k z from equation (50) and find the eigenvectors 
corresponding to each k z . Constructing general solutions for SB±, 5v±, and s from these 
eigenvectors, we find the eigenvalue condition 7(M, s' ) that must hold in order to satisfy the 
three necessary boundary conditions on the system. The boundary conditions we apply are 
Sv± = at z = ±L and the upstream boundary condition s = at z = —L. The eigenvalue 
condition thus obtained is 



for the nth order unstable mode where n — 1, 2, 3, 

A plot of constant 7 contours is displayed in Figure 3. Note that the qualitative behavior 
pointed out in the text in section 3.1 is demonstrated by this simplified model. Therefore, the 
magnetic shear is not responsible for the unexpected decrease in the stability threshold with 
shear flow; only a localization of the instability is necessary to demonstrate this characteristic. 
A quantitative comparison of Figures 2 and 3 shows that, to yield equivalent growth rates, a 
much larger entropy gradient must be supplied in the unbounded case with magnetic shear 
than in the bounded, straight-field model. We can understand the difference as follows. The 
energy required to bend the magnetic field lines slows the growth of the instability. In the 
bounded case, this energy is needed to bend the field line only within the bounded domain. 
But, in the unbounded model, the bending of the field lines occurs over a larger extent in z, 
thus requiring more energy and so more effectively suppressing the instability. 



(7 + iMk z ) [(1 - M 2 )k 2 z + 2i-fMk z + 7 2 - s' ] = 0. 



(50) 




(51) 
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3.4. Energy Conservation 

In the standard way, equations (21)-(23) can be combined, replacing 7 with the time 
derivative J^, to obtain an energy integral for the model. Converting back to 5v± and 5B± 
using equations (24) and (25), we find 



d_ 
di 



5B 2 + 5v 2 2M5B ± 5vi 



1 + z 2 



1 + z 2 



dz 



M 



s5B A 

TT~z 2 



-dz. 



(52) 



Adopting the terminology of Hayashi and Young (1987), we define the integral on the left- 
hand side of equation (52) as the wave energy of the perturbation. Note that, in the absence 
of flow (M = 0), the wave energy is constant in time. Equation (52) supplies a limit on 
the value of s' Q necessary for instability: since the 5B 2 _ and 8v\ terms are both positive 
definite, an instability can only develop for s' Q > 0. In this case, gravitational potential 
energy from the s 2 term can be harnessed to drive the kinetic energy and field line bending 
of the instability. 

By the same method as above, we find that the wave energy integral for the straight-field 
model bounded at z — +L, in section 3.3, has the form 



d_ 
di 



-L 



SBl + 5vi-*- 



dz 



-L)-5B 2 ± (L) + 



(53) 



Without flow, again, we find that the necessary condition for instability to develop is s' > 
and that the wave energy integral is constant in time. In the presence of flow, we interpret 
the terms on the right-hand side of equation (53) as follows: 5B 2 _(—L) represents the flow of 
magnetic energy into the region, —5B\(L) represents the flow of magnetic energy out of the 
region, and ^-7^ represents the flow of gravitational potential energy out of the the region. 



4. Discussion 

To study the effect of magnetic shear and shear flow on local gravitationally induced 
instabilities, we have constructed a simple model in the Boussinesq limit of ideal MHD. 
Numerical solutions to this model yield a stability diagram of the (M, s' Q ) parameter space. 
This stability diagram has three important characteristics. First, there exists a threshold 
entropy gradient for unstable mode growth, demonstrating that magnetic shear is a sta- 
bilizing influence. Second, flow serves to suppress mode growth in a substantially unstable 
regime, but near marginal stability it lessens the stabilizing effect of magnetic shear, enhanc- 
ing unstable mode growth rates and lowering the threshold entropy gradient necessary for 
instability. Third, normal modes of instability are stabilized completely as the Alfven Mach 
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number approaches one because the disturbance is swept downstream out of the unstable 
region. Analytical work corroborates these characteristics in the different regions of (M, s' Q ) 
parameter space. 

In a shear magnetic field without flow, the characteristic shape of the unstable mode 
is such that the field lines remain on the constant y' surface shown in Figure 1; hence, 
in the lab frame, the field lines must twist as they fall under gravity to remain on this 
surface. This occurs because the perturbed field line at any vertical height x must align with 
the direction of the unperturbed field at that height to facilitate interchange. Unlike the 
ordinary interchange of straight field lines, if magnetic shear is present, the field line must 
be bent to allow interchange; this is the root of the stabilizing influence of the shear magnetic 
field. Energy extracted from gravitational potential energy as the field line falls must supply 
both the kinetic energy of the moving plasma, which is frozen to the field line, and the 
energy required to bend the field line. Unstable motions are localized about the point where 
the constant y' surface is vertical (z = in Figure 1), because this point is where motions 
on that surface extract the most gravitational potential energy while minimizing field line 
bending. When plasma shear flow is included in the problem, this unstable region moves in 
the lab frame with the flow velocity Vf in the — z direction. Therefore, this characteristic 
twisting geometry propagates along z at a speed Vf = VQ l -f- determined by the local flow 
speed vq, flow shear length l v , and magnetic field shear length l B . Typical eigenfunctions for 
an entropy gradient of s' = 0.6 and Alfven Mach numbers M = and M = 0.8 are displayed 
in Figure 4. 

The modifications of the gravitational stability of a magnetized plasma due to the 
presence of magnetic shear and shear flow impact several diverse subjects. In solar physics, 
one of the key processes necessary for the success of the interface dynamo (Parker 1993; 
Charbonneau and MacGregor 1996) is the storage of an intense toroidal field in the solar 
tachocline (Spiegel and Zahn 1992) until an instability causes an isolated flux tube to rise 
into the base of the convection zone (Adam 1978; Cally 1983; Cally and Adam 1983; Hughes 
and Cattaneo 1987; Cattaneo and Hughes 1988; Cattaneo et. al. 1990; Schussler et. al. 1994; 
Caligari et. al. 1995; Gilman and Fox 1997; Barnes et. al. 1998). Shear flow is present due 
to the large differential rotation of the sun in the tachocline region (Thompson 1996; Schou 
et. al. 1998), and some magnetic shear probably exists here as well; the impact of shear 
flow and magnetic shear may alter gravitational stability within this region. There are two 
regions in galactic physics where the conditions in our model may apply: near the center 
of the galaxy, where there likely exists a region of magnetic shear in which the large-scale 
magnetic field changes from primarily azimuthal to primarily vertical (Zweibel and Heiles 
1997) and shear flow is present through differential rotation, and in the disk of the galaxy, 
where both magnetic shear and shear flow may be present in the vertical direction away 
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from the mid-plane. In examining the mechanisms for turbulence in accretion disks, Balbus 
and Hawley (1991) pointed out that a previously discovered (Velikhov 1959; Chandrasekhar 
1960, 1961) but unappreciated linear MHD instability, driven by differential rotation coupled 
with magnetic tension, would occur if the accretion disk were threaded by a weak magnetic 
field. But the Balbus-Hawley instability bends the magnetic field lines extensively; thus, it 
is stabilized by magnetic tension for large magnetic field strengths (Blaes and Balbus 1994; 
Urpin 1996; Kitchatinov and Riidiger 1997). Although the twisting interchange instability 
studied here depends not on rotation but on an entropy gradient, it has a characteristic 
geometry that minimizes field line bending and so may be important in regions of large 
field strength. As is well known within the fusion community (Roberts and Taylor 1965), 
employing a sheared toroidal magnetic field in a tokamak can help to stabilize, or at least 
suppress, ballooning instabilities (buoyant interchange instabilities driven by pressure and 
curvature forces). The idea of employing shear flow to further stabilize ballooning instabilities 
has gained much attention in the past decade (Waelbroeck and Chen 1991; Hassam 1991, 
1996, 1999; Miller et. al. 1995). Our work demonstrates that, near marginal stability, 
shear flow actually lessens the stabilizing effect of magnetic shear, lowering the threshold 
entropy gradient required for instability and enhancing the unstable mode growth rate. Full 
stabilization of the plasma will occur only if the Alfven Mach number of the plasma flow 
(as defined in our transformed coordinates) exceeds one. The local nature of the instability 
examined in this work means that our treatment may apply locally in more complicated 
environments. The instability may behave as a traveling "wave packet" which moves with 
the intersection of shear magnetic field lines (where the constant y' surface is vertical); the 
disturbance may move in and out of unstable regions, with the perturbation growing where 
conditions are unstable and decaying in stable regions. 

We have extended the model to include compressibility. In this case, the behavior 
is governed by a system of five coupled, first-order ordinary differential equations. Five 
parameters are necessary to describe the system: the Alfven Mach number, the plasma /3, the 
density gradient, the pressure gradient, and the magnetic field gradient. Over some portions 
of this five-dimensional parameter space, the growth rate eigenvalue 7 does indeed become 
complex. Producing a simple answer from this more complicated model is quite difficult. Our 
current research is addressing this difficulty But, the magnetic-buoyancy instability (Parker 
1979), which depends on compressibility, cannot be investigated without employing this more 
detailed treatment. We will then be able to relate our work quantitatively to applications 
such as the stability of magnetic fields in the solar tachocline. One final possible extension 
of this research is an investigation of the nonlinear behavior using an initial-value code. 

The research was performed under appointment of Greg Howes to the Fusion Energy 
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A. Twisting Modes vs. Fourier Modes 

In the sheared magnetic field, the local unstable modes can be written in two ways 
(Roberts and Taylor 1965); we call them twisting and Fourier modes. For the reader's 
convenience, we summarize the essence of the argument here. In the twisted field line 
coordinates, we obtain the twisted mode where a perturbed quantity $(x, y, z, t) is given by 

$ T = $ T (z, x)e iky ' + * = ® T (z, x)e ky ~ lk ^ +lt (Al) 

with $^ localized in z and varying weakly in x compared to k~ l . But since the origin in z 
is arbitrary, we can also write 

^ = $ T ( Z _ Zo , x)e lky ~ lkXA ^ 1+lt . (A2) 

Thus, there are an infinite number of twisting modes, each with a different origin of the 
twist. We can construct a mode that does not depend on z — a Fourier mode — by integrating 
equation (A2) over z 

& T (z - zo, x, y, t)dz = $ F I—, x) e^ 1 , (A3) 



where 

* F = f <5 T (z',x)e-^ z ' dz '. (A4) 

Thus, the Fourier modes and the twisting modes are related by a Fourier transform. Clearly, 
the Fourier mode can be made by "adding" twisting modes together (equation [A3]), or 
vice versa, using the Fourier inversion theorem on equation (A4); see Cowley et. al. (1991) 
for pictures of this superposition. The Fourier modes are narrowly localized in x — typically 
Ax ~ where Ax is the x width of the Fourier mode and Az is the z width of the twisting 
mode. Note this x localization of the Fourier mode is narrow compared to the x variation 
of the twisting mode. In this paper, we have taken the twisting mode representation for two 
reasons: first, the role of the flow, we believe, is more intuitive in this picture; and, second, 
the twisting modes are finite in z extent and therefore represent more easily the evolution 
of an initial value problem. 
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B. Asymptotic Analysis as M — > 1 

To demonstrate stabilization as the Alfven Mach number approaches one, as seen of 
region (II) of Figure 2, we perform an asymptotic analysis of our model in the limit M — > 1. 
For a small dimensionless parameter e, we quantify the order of this limit as 1 — M ~ 
0(e 2 ). We expect the instability growth rate to be 7 ~ O(e). Identifying the terms in the 
dimensionless system of equations (21)-(23) for reference, we have 

-(1-^)^ = -7^ + +(1 + M)^4_ 

(1) (2) (3) (4) 1 ' 



(5) (6) (7) (8) 

M£ = -is - s -ki^m ^ 



;b2) 



dz 1° 2 (l+z 2 )V2 "T 2 (l+z 2 )V2 ■ mg\ 

(9) (10) (11) (12) 

The boundary conditions demand A + ,A_,s — > as |z| — > 00. Asymptotic solutions in the 
M — > 1 limit can be found in the four regions along z displayed in Figure 5. Below we find the 
solutions for each of these regions and, by matching the solutions between these regions, we 
obtain an eigenvalue condition on the growth rate demonstrating stabilization as the Alfven 
Mach number approaches one. First, we obtain two reductions of equations (B1)-(B3); one 
over a boundary layer where \z\ <C e -1 , and the other over an outer region where \z\ >> 1 . 
Then, we present the solutions in each of the four regions in Figure 5. 

As the Alfven Mach number approaches one, regions (1), (2), and (3) of Figure 5 behave 
like a boundary layer region: we expect derivatives to be large and thus take ^ ~ 0(e^ 1 ). 
We treat z ~ 0{1) over these three regions. Balance of the dominant terms (9) and (11) 
in equation (B3) shows that s ~ 0(eA + ). In turn, this demands that terms (5) and (8) 
must balance in equation (B2), yielding the ordering A_ ~ 0(e 2 A + ). Adopting the specified 
ordering allows us to drop terms (6), (7), (10), and (12) from equations (B1)-(B3); term (3), 
although one order in e smaller than the other terms in equation (Bl), will contribute in the 
regions for z > 1, so we retain it in order to be certain that our first order correction in region 
(3) is valid. The remaining equations can be combined to a single third-order equation in 
A_ and simplified by the substitution z = sinh 9 to obtain 

2M(l-M)d 3 A„ 2 7 M , d 2 A^ dA_ 

In region (4), we find the smoothly varying outer solution over which z ~ C^e" 1 ). We 
expect that j- z ~ 0(e) here. By insisting that terms (9) and (10) balance with term (11) 
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in equation (B3), we find that s ~ 0(A + ). Similarly, by balancing terms (5) and (6) with 
term (8) in equation (B2), we obtain ~ 0(A + ). This ordering allows us to drop terms 
(1) and (7) in equations (B1)-(B3). Approximating (1 + z 2 ) 1 ^ 2 ~ z, substituting A_ = zB_, 
and combining the equations, we obtain the second-order equation 



M(l + M)z 2 ^ r +[(l + 2M) 7 z 2 +3M(l + M)2-^ 



dB. 
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dz 

+ [ 7 V + (1 + 3M) 7 2 + M{1 + M) — s' ] B_ = 0. (B5) 

To find the solution in region (1), we assume an eikonal solution for equation (B4). 
Neglecting the trivial constant solution, we find two independent solutions of the form 

A - " T7? ex P ( tu, 1 T7\ [ z + ( 1 + z2 ) 1/2 } ) (B6) 

(l + ^)3/4[( 1 + z 2)l/2 + z ]V2 ^\2(1-M) 1 V 

A - ~ Z \ilVyT exp (w^m) [z ~ {l + z2f,2] + ^ " (1 + ^ 2)1/2 • (B7) 

To get the behavior for \z\ 3> 1, we can expand (1 + z 2 ) 1 / 2 ~ \z\ + l/2|z|. For region (1), we 
note that z < and, retaining only the dominant terms, we obtain solutions of the form 

(B8 > 

1 - ( ->> • (B9) 



The boundary conditions impose that A_ — > as z — > — oo, so our solution in region (1) 
must be entirely of the form of equation (B9), a growing solution in the +z direction. To 
determine the behavior of equations (B6) and (B7) in the overlap with region (2), we take 
the limit \z\ <C 1 and approximate (1 + z 2 ) 1 ^ 2 ~ 1 + z 2 /2 to get the two solutions 

A -~^{w^) [z + z * l2] ) (B10) 
- 4 -~ exp (^?) [; "" /2] )- (B11) 

One of these solutions must smoothly match onto the solution for region (2). 

Region (l)'s solution will be valid as we move in the +z direction until the eikonal 
approximation, ^> 1, breaks down. The failure of this condition occurs in region (2) of 
Figure 5 where z ~ 0(e). For region (2), we expand equation (B4) about z — 0. For z C 1, 
cosh# ~ 1 and tanh# ~ 9 where ^ < 1, so we can drop the last term of equation (B4). 
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Letting / = and using the integrating factor f — B_ exp ( 2 T (i-m) ) ^° simplify the result, 
we obtain the equation 



d 2 B_ r 7 2 cosh 2 g 7sinhg s' 1 

d# 2 U(!- M ) 2 2(1 — M) 2M(1-M)J 1 ' 

Neglecting the central term in the coefficient of B_ because it is an order e smaller than 
the other terms, we can cast equation (B12) in the form of Hermite's equation for which 
the solutions are well known. Therefore, for the n = Hermite polynomial, the solution in 
region (2) is 

A- - few (^Tm) V ~ z ' 2 / 2 )) dz 'i (B13) 
and the eigenvalue condition on the growth rate imposed by Hermite's equation is 

7 2 = 2S> ° {1 ~ M) - 2(2n + 1)7(1 - M) (B14) 

for the nth Hermite polynomial. Thus, we find that the solution in region (1) given by 
equation (Bll) matches smoothly onto our solution in region (2). 

In region (3), we once again assume an eikonal solution for equation (B4) and find the 
two solutions given by equations (B6) and (B7). To match with region (2), find the \z\ <C 1 
limit of these equations, yielding once more equations (BIO) and (Bll); we observe that 
equation (Bll) for the small z limit of region (3) matches solution (B13) for region (2). In 
the \z\ 1 limit of equations (B6) and (B7) for region (3), we obtain the solutions (for 
z > 0) 

A_~l e xp(-^-) (B15) 



z 2 L \l-M, 

A -~ exp (i(i^)i)' < B16 > 

To continue our asymptotic solution, we must smoothly match one of these region (3) solu- 
tions to the solution for region (4) in the overlap around z ~ 0(e~ 1 / 2 ). 

In region (4), we assume eikonal solutions for equation (B5) in the limit z — > 0. The 
two solutions found are 

A_ ~ z (B17) 

A -~^(^m)- < B18 > 

Hence, we can match the solution given by equation (B16) in region (3) with the solution 
given by equation (B18) in region (4) if 

S '° (B19) 



4(1 — M) 2 7 M' 
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But, this is identical to the lowest order of the eigenvalue condition, equation (B14). To 
complete our asymptotic solution, we must find that a solution to equation (B5) in the limit 
z — > oo which satisfies the boundary condition that A_ — > as z — > oo. In this limit, the 
two solutions take the form 



Both of these solutions satisfy the boundary condition as z — > oo. 

Now that we have seen that it is possible to construct a complete asymptotic solution 
in the limit M — > 1, let us examine this solution more closely. Beginning in region (1) at 
the left of Figure 5, the boundary conditions demand that the solution must be solely of 
the form of equation (B9). But, as behavior in regions (1), (2), and (3) is governed by 
equation (B4), the eikonal approximation must break down in region (2) in order for the 
solution in region (1) to convert to the solution given by equation (B16) in region (3) so 
that smooth matching may be accomplished with solution (B18) in region (4). The failure of 
the eikonal approximation around z = yields a reduction of equation (B4) to the Hermite- 
type equation (B12). The requirement that a solution to this equation exist imposes the 
eigenvalue condition, equation (B14). This single condition can also be used to smoothly 
match solution (B16) in region (3) to solution (B18) in region (4). Finally, region (4) is 
governed by equation (B5). In the limit z — > 0, this equation yields a matching solution in 
the overlap with region (3); and, in the limit z — > oo, it provides two solutions which both 
satisfy the boundary conditions as z — > oo. Therefore, the single condition necessary to 
find a smooth solution which satisfies the boundary conditions is the eigenvalue condition, 
equation (B14). To lowest order, this condition can be written in a more recognizable form 
in the limit M — > 1 as 



Thus, in region (II) of Figure 2, where 7 ^> 1—M, the behavior is clearly demonstrated — that 
stabilization occurs as the Alfven Mach number approaches one. 
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(B21) 



7 2 ~ Sq(1 — M 2 ). 
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Fig. 1. — The geometry of the shear magnetic field as well as the applied shear flow (shown 
as vo) is shown. The constant y' surface is represented by the magnetic field lines (solid 
lines) and the dashed lines. 
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Fig. 2. — Stability diagram for the Boussinesq limit: Contours of constant normalized growth 
rate 7 are plotted over the space of entropy gradient s' vs. Alfven Mach number M. The 
stable parameter regime is denoted by hashing. The diagonal dotted line denotes 7 = 1 — M, 
separating unstable region (I), where flow enhances instability growth, from unstable region 
(II), where flow suppresses the instability. 
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Fig. 3. — The stability diagram for the straight-field case as described by equation (51) in 
section 3.3. Note that the qualitative features of the diagram are similar to those in figure 2. 
A value of L = n was chosen to plot this diagram. 
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Fig. 4. — Eigenf unctions for entropy gradient s' = 0.6 and M = (left plot) and M = 
0.8 (right plot). The three functions are A + (solid line), A_ (dotted line), and s (dashed 
line). Note that for the M = 0.8 case, where the plasma flow is in the +z direction, the 
eigenfunctions grow quickly with steep gradients and diminish slowly as you move from left 
to right. 
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Fig. 5. — Regions denned for the asymptotic solution of A_ in the limit M — > 1. 



